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Abstract 

We study the dynamics of solitons in a Landau-Lifshitz equation describing the mag- 
netization of a three-dimensional ferromagnet with an easy axis anisotropy. We 
numerically compute the energy dispersion relation and the structure of moving soli- 
tons, using a constrained minimization algorithm. We compare the results with those 
obtained using an approximate form for the moving soliton, valid in the small mo- 
mentum limit. We also study the interaction and scattering of two solitons, through 
a numerical simulation of the (3+l)-dimensional equations of motion. We find that 
the force between two solitons can be either attractive or repulsive depending on 
their relative internal phase and that in a collision two solitons can form an unstable 
oscillating loop of magnons. 
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1 Introduction 



In the continuum approximation the state of a ferromagnet is described by a three 
component unit vector, 0(x, t), which gives the local orientation of the magnetization. The 
dynamics of the ferromagnet, in the absence of dissipation, is goverened by the Landau- 
Lifshitz equation 

a=^ X ^ (L1) 

where E is the magnetic crystal energy of the ferromagnet. We have chosen units in which 
the spin stiffness and magnetic moment density of the ferromagnet are set to one. 

The case we study in this paper is that of a three-dimensional ferromagnet with isotropic 
exchange interactions and an easy-axis anisotropy, in which case the energy is given by 

d i( f> + A{1 - 4>D) rf 3 x (1.2) 

where A > is the anisotropy parameter and we choose the ground state to be <fi = 
(0i,0 2 ,0 3 ) = (O,O,l) = e 3 . 

In this case the Landau-Lifshitz equation becomes 

^ = <£x(W + ^ 3 e 3 ). (1.3) 

This equation has finite energy, stable, exponentially localized solutions known as mag- 
netic solitons ||. In section 2 we review the properties of stationary magnetic solitons and 
in section 3 we numerically compute solutions describing moving solitons, display their 
energy dispersion relation and discuss their structure. Finally, in section 4, we perform 
numerical simulations of the time dependent equations of motion to investigate the interac- 
tion and scattering of two solitons. We find that the force between two solitons depends on 
their relative internal phase, and that during a collision two solitons can form an unstable 
loop of magnons which subsequently decays into solitons. 



2 Stationary Solitons 



In addition to the energy flO), the Landau-Lifshitz equation ( |1.3|) has two other con- 
served quantities. These are the number of spin reversals, N, and the momentum P, given 
by! 

N= /(l-0 3 ) rf 3 x (2.1) 



and 



1 r 

Pi = -fijk J Xjtkim4> ■ di4> x d m <p d x. 



(2.2) 



In the quantum description, N counts the number of quasi-particles in the magnet, that 
is, it may be interpreted as the magnon number. 
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Figure 1: The number of spin reversals, N, as a function of the precesion frequency u. The 
solid curve is the numerical result and the dashed curve is the approximate formula 



In this section we consider only stationary solitons, ie P = 0, whose properties have 
been discussed in ref. [H]. There are no static solitons, but stationary solitons have a time 
dependence which resides only in the constant motion of an internal phase. Explicitly, the 
stationary soliton has the form 

= T^7l( 2 / cos M)' -2/sinM), 1 - f) (2.3) 
1 + f 2 

where u is the frequency and f(r) is a real profile function with the boundary conditions 
/(oo) = 0, /'(0) = 0. The resulting equation for the profile function is 



r = _?f + vf + Afii-r) _ uf ( , 4) 

r 1 + / 



For large r this equation may be linearized to give the asymptotic solution 



/ ~ B exp(-rvM - uo) (2.5) 

from which it can be seen that a soliton solution exists only for u < A. It can also be shown 
that uj > and hence an anisotropy term is vital for the existence of these stationary 
solitons. From now on we set A = 1 so that < u < 1. There is a one-parameter 
family of stationary solitons parameterized by either the frequency uo or the number of 
spin reversals N. In figure 1 we plot the relation between these two quantities, obtained 
by solving equation ( |2.4|) using a shooting method (solid curve), and in figure 2 we display 
(solid curve) the energy per spin reversal of the soliton, E/N as a function of N. 
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Figure 2: The energy per spin reversal, E/N, as a function of the spin reversal N. The 
solid curve is the numerical result and the dashed curve is the approximate formula (|2.6|) . 



In the limit of small u (which corresponds to large N) the solution is approximately 
given by / = e R ~ r , where R ^> 1 is a core radius which marks the boundary between 
the two vacua <f> = ±e 3 inside and outside the soliton. In this case the soliton may be 
thought of as a spherical domain wall. Using this simple approximate solution the following 
asymptotic formula may be derived, which are valid for large N 

N=% E= 3 -^ = (72nN^. (2.6) 

The dashed curves in figures 1 and 2 are obtained from these asymptotic formulae and fit 
the numerical results well. 



3 Moving Solitons 

To compute moving solitons is a more difficult task than in the stationary case and 
can not be reduced to simply solving an ordinary differential equation. There is an ansatz 
which is consistent with the equations of motion and describes a soliton which moves at 
constant velocity v and rotates in internal space with frequency u. Explicitly the ansatz 
reads 

0i+# 2 = (0i(x- vt) +*0 2 (x-vt))e-^, 3 = 3 (x-vt). (3.1) 

Substituting this ansatz into the equation of motion ( |1.3j ) leads to a partial differential 
equation for <fi which is not compatable with a spherically symmetric solution for non- 
zero v. This partial differential equation has a variational formulation, which is useful 
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in computing its solutions. Let E(P,N) be the minimal value of E for fixed values of 
the momentum P and number of spin reversals N. Then the solution of this constrained 
minimization problem is precisely the function (j> corresponding to the values H 



dE 
dN 



dE 



dp 



N 



(3.2) 



We compute moving soliton solutions by numerically solving this constrained minimiza- 
tion problem. A similar computation has been performed in the case of an isotropic magnet 
(A = 0) 0, but the results in this case are qualitatively different. For example, we have 
already seen that without anisotropy there are no stationary solutions, and in fact moving 
solitons exist in the isotropic case only for sufficiently large momenta ||. 

A moving soliton has an axial symmetry in the plane perpendicular to its momentum 
P. Since the Landau-Lifshitz equation is 5*0(3) invariant we can, without loss of gener- 
ality, choose the momentum to be in the x% direction, P = (0,0, P), so that the soliton 
moves along the x% axis and has an axial symmetry in the X\X2 plane. We use cylindrical 
coordinates, with p - 



I x\ + a?2, and discretize the energy using second order differences 
and a grid of size 50 x 100 in the p, X3 plane. The energy is minimized using a simulated 
annealing algorithm || |J and the constraints on iV and P are imposed using a penalty 
function method. 




P/N 



Figure 3: The energy, E, as a function of the scaled momentum, P/N 2 / 3 , for the three 
values N = 204, 317, 542. The crosses represent the numerical data and the curves are 
obtained from the small momentum approximation. 



In figure 3 we plot (crosses) the energy dispersion relation, E(P, N) obtained from these 
computations for the three values N = 204, 317, 542. The energy increases with N so the 
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Figure 4: The spin reversal density, 1 — 03, in the x 2 = plane, for the soliton with 
iV = 317, and four values of the scaled momentum, P/N 2 ^ 3 = 0, 6.8, 13.6, 68.2 



data need not be labelled individually. Note that the rather unobvious choice for these 
three values is due to the fact that in the zero momentum case they correspond to the 
frequencies u = 0.7, 0.6, 0.5. The soliton velocity, if required, can be read-off from this data 
by making use of formula ( |3.2| ). 

We now describe an approximate initial condition for the field </> to describe a soliton 
with any given values of N and P and assess its accuracy by computing its energy dispersion 
relation and comparing with the results displayed in figure 3. 

It was noted in ref.JTJ that in the two-dimensional case a moving soliton could be 
generated from initial conditions which consist of performing a space dependent phase 
transformation on the stationary soliton solution, although no relation was given between 
the parameter which appears in this transformation and the momentum of the resulting 
soliton. In the following we study the same type of transformation in the three-dimensional 
case, although we have a simple explicit form for the parameters N, P of the corresponding 
solution and this will allow us to compute the energy dispersion relation to compare with 
the previous numerical results. 

Our approximate initial conditions are given by 

^ = T ^(2/cos(Px 3 /iV),2/sin(Px3/iV), 1 - f) (3.3) 

where /(r) is the stationary soliton profile function with N spin reversals, which we have 
already introduced earlier. Note that since 03 is the same as in the stationary case then the 
number of spin reversals of this field is indeed N. We now show that the momentum of this 
field is P = (0, 0, P). By symmetry it is trivial to show that Pi = P 2 = 0. Substituting (|3.3j ) 
into the momentum equation ( |2.2|) , using spherical variables and performing the angular 
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part of the integration we have 



16vrP r°° ff'r 3 , 8vrP r°° f 2 r 2 1. /V 3 , ,„ 

P3 = -^v-io (Tttt dr= ^rL ^n'^^n ] ( ) 

where we have used the boundary conditions on f{r) to set the total derivative term to 
zero and recognized in the remaining term the expression for N. 

Having verified that this configuration has N spin reversals and momentum P we can 
now compute its energy and derive its dispersion relation. Substituting the field (|3.3|) into 
the expression for the energy (|1.2|) and performing the angular integrals we find 



E = E {N) + P 2 E l (N), (3.5) 

where 

roo (fa , f2\ 2 g ,oo f2 2 



The dispersion relation obtained from equation (|3.5| ) is also displayed in figure 3 (solid 
curves) for comparison with the full numerical results. In ref . it is noted that a magnon 
loop occurs for P iV 2//3 , and from figure 3 it can be seen that the ansatz is only a good 
approximation if the momentum stays below this range, since otherwise it begins to lose 
its accuracy and has considerably more energy than the true solution. 

In figure 4 we plot, in the X2 = plane, the spin reversal density, 1 — 03, for the soliton 
with N = 317 and four different values of the scaled momentum, P/N 2 / 3 = 0, 6.8, 13.6, 68.2 



Recall that in the spherical ansatz (3.3) the spin reversal density is assumed to be 



the same as in the stationary case. Therefore by examining how the spin reversal density 
changes with increasing momentum we can assess the validity of the ansatz. By a com- 
parison of figures 4a and 4b, it can be seen that the spherical assumption in the ansatz 
is substantially violated for P ~ N 2 ^ 3 , with the true soliton being stretched in the plane 
orthogonal to its motion. For very large momenta, P N 2 ^ 3 , there is a transition from 
a single lump of magnons into a magnon loop. Already at P/N 2 ' 3 = 13.6, fi gure 4c, one 
can see that the soliton is sufficiently distorted so that the maximum of the spin reversal 
density is no longer on the X3 axis. For P/N 2 ' 3 = 68.2, fi gure 4d, the soliton has now 
clearly formed two distinct components, which when we recall the axial symmetry in the 
X1X2 plane, means that the spin reversal density is now localized around a circle in a plane 
perpendicular to its motion. In the isotropic case all the solitons have this structure 
and are known as magnetic vortex rings. 

In summary, we have seen that low momentum solitons can be described by a spherical 
ansatz but as the momentum increases there is a transition from lump-like solitons to 
ring-like solitons. 
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4 Multi-Soliton Interactions 



In this section we discuss the results of a numerical evolution of the full time- dependent 
Landau-Lifshitz equation (|1.3|) in order to investigate the interaction and scattering of two 
solitons. 

Two initially stationary and well-separated solitons have an axial symmetry about the 
line connecting them, and we make use of this symmetry in our numerical evolution code. If 
we take two solitons on the X3-axis, each with momentum only in the X3 direction, then both 
the initial conditions and the equations of motion have an axial symmetry in the X\Xi plane. 
We could use cylindrical coordinates and evolve the equations of motion in these variables 
but this type of implementation can often lead to numerical instabilities associated with 
coordinate singularities, particularly along the a^-axis. We therefore employ the 'Cartoon' 
method recently proposed in ref.[ffl] in the context of solving the axisymmetric Einstein 
equations. In this approach the equations are evolved in cartesian coordinates, but on a 
thin slab, that is the grid has size M x 3 x M so that only 3 points are used in the X2 
direction. In the centre of the slab, which is the plane £2 — 0, the equations of motion are 
evolved in a standard manner; we employ a fourth-order Runge-Kutta method for the time 
evolution and space derivatives are approximated by second order finite differences. Away 
from this central slice the field values are determined by making use of the axial symmetry 
to map to an equivalent point in the central slice, x<i = 0. Since the field values in this 
central slice are only known on a discrete lattice a one-dimensional interpolation algorithm 
must be employed, which we take to be a simple second order Lagrange interpolant. The 
results in this section were obtained using grids of size M = 201, with time and space steps 
given by At = 0.025, Ax = 0.3. 

To obtain an initial condition consisting of two well-separated solitons we use the fol- 
lowing nonlinear superposition rule. Given the S^-valued field 4> corresponding to the first 
soliton we construct the associated Riemann sphere variable obtained by stereographic 
projection of the S 2 - valued field onto the complex plane. We construct a similar Riemann 
sphere field for the second soliton and obtain the combined S 2 field by addition of the two 
Riemann sphere fields followed by an inverse stereographic projection. Note that before the 
two fields are combinded a relative phase, a, can be introduced corresponding to rotating 
the 0i, <p2 components of the first soliton through an angle a. As we shall see, this relative 
phase has a marked influence on the force between the two solitons. 

As an initial condition we take two solitons, each corresponding to a stationary soliton 
with N = 317, placed on the x 3 -axis at the positions x 3 = ±7. Figure 5 diplays the 
subsequent evolution of their relative separation for the cases when there is no relative 
phase, that is a = 0, (bottom curve) and when the two solitons are exactly out of phase, 
that is a = 71, (top curve). In computing the separation of the two solitons the position 
of each soliton is defined as the point at which the spin- reversal density, 1 — 3 , takes its 
maximum value, and this is calculated using a quadratic interpolation around the lattice 
site at which it takes its maximal value. 

From figure 5 we see that when there is no relative phase the two solitons attract and 
form a single larger soliton, which is perhaps what is expected from the fact that the energy 
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Figure 5: The time evolution of the separation between two initially stationary solitons. 
The relative internal phase, a, is set to (bottom curve) and 7r (top curve). 



per spin reversal is a decreasing function of N as shown in figure 2. The slight peak before 
the solitons merge in figure 5 is a result of the highly nonlinear deformation the solitons 
suffer as they combine. In contrast, when the solitons are out of phase they repel and 
quickly settle to a state in which the two solitons are moving away from each other at a 
constant speed. Thus we have demonstrated that the force between two solitons can be 
attractive or repulsive depending on their relative phase. 

If the relative phase is something other than or 7r then the situation is even more 
complicated. In this case the ultimate fate of the solitons, that is whether the final configu- 
ration consists of one or two solitons, depends on the relative phase and initial separation of 
the solitons. However, the picture is not as simple as a mere attraction or repulsion, as the 
solitons can initially move towards each other but then turn around and ultimately repel. 
Furthermore, during this interaction there can be a magnon exchange, so that the number 
of spin reversals of one soliton can increase while that of the other decreases. Qualitatively 
similar results have been found and studied in great detail H for solitons known as Q-balls, 
where the role of the number of spin reversals is played by a Noether charge. Q-balls occur 
in relativistic field theories, in which the dynamics is second order in time, rather than 
the first order dynamics of the Landau-Lifshitz equation, so it may seem a little surprising 
that there should be similarities. However, Q-balls have a time-dependent internal phase, 
very similar to the precession of the magnetization for the magnetic solitons studied in 
this paper, and as argued in ref.0, to which we direct the reader for further details, the 
novel dynamics can be mainly attributed to this fact. As we discuss below, there are other 
common features between the dynamics of magnetic solitons and Q-balls. 

For two solitons which individually have no momentum and are in-phase we have seen 
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that the two solitons merge and, after some oscillations, settle down to a single soliton 
which has a value of N which is approximately the sum of its constituents. A small 
amount of energy is radiated during these oscillations, in the form of magnons which are 
released to infinity, and can be dealt with numerically by applying absorbing boundary 
conditions at the edge of the grid. As we now demonstrate, if the solitons initially have 
some momentum then a much more complicated evolution takes place. 



j -> O ,:j 

1 2 3 
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Figure 6: Isosurface plots of 1 — 3 = 1.5 at times t = 0,20,30,65,80,115 during the 
collision of two solitons. 

In figure 6 we display a fully three-dimensional isosurface plot corresponding to the 
surface where 1 — 03 = 1.5, at six different times t = 0,20,30,65,80,115. The initial 
conditions were created from two N = 317 solitons, with zero relative phase, placed at the 
positions 23 = ±7 and with momentum P = =f158. The individual moving solitons were 
constructed using the ansatz ( |3.3|) , which is an adequate approximation at this momentum. 
We see from figure 6 that the two solitons collide and form a single lump (figure 6.2) which 
then expands to form an axially symmetric loop of magnons, that is, a toroidal region 
where the spin-reversal density is concentrated (figure 6.3). This magnon loop initially 
expands until it reaches a critical radius, which increases with the initial momentum of the 
solitons, and then contracts to produce two solitons which move apart along the initial line 
of approach (figure 6.4). The solitons then attract once more, reforming the loop (figure 
6.5), which again collapses, with the solitons again attempting to separate (figure 6.6). 
This process continues through several cycles, with the loop forming at a slightly smaller 
radius each time, until eventually the configuration settles to a single large soliton with 
no momentum. Note that the loop which forms in this process is not the same kind of 
structure as the moving vortex ring, since it has no momentum perpendicular to the plane 
of the loop. If the initial momentum of the solitons is increased slightly (for example, 
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with P = 317, and all other initial conditions the same) then the first loop formed is of 
sufficient radius that its collapse results in the two solitons which emerge from its decay 
having sufficient momentum that they never recombine and instead travel out to infinity. 

This phenomenon is the three-dimensional realization of the two-dimensional process 
described in ref. |7j where two solitons scatter at right angles to the initial line of approach. 
In the three-dimensional case the axial symmetry prevents the two solitons from scattering 
at right angles and results instead in the formation of a magnon loop. 

This kind of loop formation also appears in relativistic Q-ball dynamics and the quali- 
tative features are very similar ||. 

In summary, we have seen that the interaction of multi-solitons, even at relatively 
low momenta, is a very complicated process with forces which depend on relative internal 
phases, and novel features such as the production of unstable loops. Clearly these processes 
require further study and hopefully a more analytical understanding will emerge. 

5 Conclusion 

We have used several different numerical techniques to study the dynamics and inter- 
action of magnetic solitons in a three-dimensional ferromagnet with easy-axis anisotropy. 
We have computed moving solitons using a minimization algorithm and compared the re- 
sults to those of a simple radial ansatz, which we have shown is a good approximation for 
low momenta. However, for large momenta there is a transition from lump-like solitons to 
ring-like solitons, where obviously the radial ansatz fails badly. We have found that the 
interaction between two solitons has a strong dependence on their relative phase and that 
the collision of solitons can be highly non-trivial and lead to the formation of unstable 
magnon loops. Finally, we have observed that many of the features found for magnetic 
solitons are qualitatively very similar to those of Q-balls, despite the fact that the dynamics 
of the latter is a second order in time relativistic system. 
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